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Abstract 



We analyze with the tools of lobe dynamics the velocity field from a numerical 
simulation of the surface circulation in the Northwestern Mediterranean Sea. We iden- 
tify relevant hyperbolic trajectories and their manifolds, and show that the transport 
mechanism known as the 'turnstile', previously identified in abstract dynamical sys- 
tems and simplified model flows, is also at work in this complex and rather realistic 
ocean flow. In addition nonlinear dynamics techniques are shown to be powerful 
enough to identify the key geometric structures in this part of the Mediterranean. In 
particular the North Balearic Front, the westernmost part of the transition zone be- 
tween saltier and fresher waters in the Western Mediterranean is interpreted in terms 
of the presence of a semipermanent "Lagrangian barrier" across which little transport 
occurs. Our construction also reveals the routes along which this transport happens. 
Topological changes in that picture, associated with the crossing by eddies and that 
may be interpreted as the breakdown of the front, are also observed during the simu- 
lation. 
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1 Introduction 



Ocean water masses of different origins have distinct contents of salt, heat, nutrients and 
chemicals. Currents transport them, and energetic mesoscale features are responsible for 
most of their mixing with surrounding waters. 

Vortices are the most well studied of such structures. Frequently they are long lived, 
and their cores remain relatively protected from the neighboring areas, so that water 
trapped inside could maintain its biogeochemical properties for long time, being trans- 
ported with the vortex. In steady horizontal velocity fields, the centers of vortices are 
readily identified as elliptic points in the streamfunction level sets. The presence of closed 
streamlines around them is the mathematical reason for the isolation of the vertex core 
from the exterior fluid. 

When the velocity field changes in time, closed streamlines are replaced by more 
complex structures, some of which can be related in idealized cases to the Kolmogorov- 
Arnold-Moser tori of dynamical systems theory. For slowly varying velocity fields vortex 
cores remain coherent during some time, but there is vigourous stirring of the surround- 
ing fluid that finally leads to water mixing. In order to understand this mixing process, 



one should focus on features of the velocity field diffe rent from the e 



terizing the vortex cores. Since some time ago (|Tu et al 



2003 



Ide et al 



iptic points charac- 



2002; Mancho et al 



2004j) hyperbolic trajectories (trajectories with saddle-like stability properties which are 



solutions to a dynamical system) have been recognized as the structures responsible for 
most of the stretching and generation of intertwined small scales that finally leads to mix- 
ing. In particular there are Distinguished Hyperbolic Trajectories (DHTs), characterized 
by its special persistence as compared with the other hyperbolic structures, that act as the 
organizing centers of the fluid stirring processes. 

Despite the great amount of attention devoted to the identi fication and characteriza- 



tion of vortices and their dynamics in ocea nographic contexts (jOlson 



2002; 



Ruiz et al. 



2002; 



Isern-Fontanet et al. 



1991 



Puillat et al 



2006), there are few studies focussed on hy- 



perbolic objects, and most of them in idealized settings. A possible reason for this may be 
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that the intrinsic instability of trajectories close to hyperbolic points makes more difficult 
their identification, in contrast with the recurrent character of trajectories in vortices that 
allow their tracking for long times from in situ and satellite measurements. In addition, 
whereas many aspects of vortex dynamics can be analyzed in a Eulerian framework, hy- 



perbolic trajectories are defined by their 



can they be fully characterized (|lde et al 



agrangian characteristics and only in this frame 



2002; 



Manchoetal 



2004). 



In this Paper we identify relevant hyperbolic trajectories from the surface velocity field 
of the Western Mediterranean Sea, obtained from a three dimensional model simulation 
under climatological atmospheric forcing. Our aim is to show that transport mechanisms, 
in particular the so called / turnstile / mechanism, preyiously identified in abstract dynam- 
ical s ystems (jChannon and Lebowitzl Il980t Bartlettl Il982l MacKav et al.l ll984L Iw iggins 



1992 ), and discussed in the context of rather simple model f lows (|Rom-Kedar et al 



Beigie et al.lll994 : 



Samelson 



1992; 



199 



0; 



Duan and Wiggins! Il99q) , are also at work in this com- 



plex and rather realistic ocean flow. More broadly, nonlinear dynamics techniques are 
shown to be powerful enough to identify the key geometric structures in this part of the 
Mediterranean. 

Western Mediterranean surface layers (up to a d epth of about 150 m) contain Modified 
Atlantic Waters of rather different characteristics iMilloJ 1999). Fresher waters (salin- 
ity about 36.5 psu) recently entered from the Atlantic occupy the Algerian Basin at the 
South and older and saltier waters (salinity above 38 psu) occupy the Northern part of 
the area. The dynamics of the contact zone between these two surface water masses, 
and eventually their mixing, is important to understand the physical and biogeochemical 
properties of the Western Mediterranean. We focus on one of the oceanographic struc- 
tures known to be of importance in these processes, the so-called North Balearic Front 



^Lopez- 



Garria et al. 



1994; 



Pinot et al. 



lese p 
19951: 



Millot 



1999). It extends roughly along the 



Southwest-Northeast direction at the North of the Balearic Islands, with significant dis- 
placements and deformations. It is characterized by a strong salinity jump of about 0.6 
psu (37.4 psu to the South and 38.0 to the North) down to 150 m depth. This identifies the 
front as the main transition zone between the two water masses. In Winter a weak but 
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detectable temperature gradient of 0.5-1 K/5 km can be observed in satellite images. 

After selecting an interval of time in our simulation during which the front is well 
formed, we explore the transport properties of the surface velocity field in the region, 
and find the relevant Lagrangian structures, hyperbolic points and their manifolds, re- 
sponsible for the establishment and permanence of the front. The location of the front is 
identified as a "Lagrangian barrier" across which transport is small (as quantified with 
the tools of lobe dynamics), and occurs via filaments that entrain water along transport 
routes that we identify. The presence of eddies strongly affect the Lagrangian structures 
in a process that can be interpreted as the break-down of the front. 

The Paper is organized as follows: In Section |2] we introduce some of the dynamical 
systems concepts that will be used in the following. Section |3] describes our numerical 
ocean model and Section |4] addresses the adaptation of the standard algorithms to the 
kind of oceanographic data provided by the model. Section |5] contain our main results 
and the conclusions are summarized in Section |6j 



2 Distinguished Hyperbolic Trajectories and their Mani- 
folds 

In recent years there have been many applications of the dynamical systems approach to 
transport in oceanographic flows. R ecent reviews are Tones and Winkler (|2002|) : Iwigginsl 
(2005); ISamelson and W iggins (2006). In this section we describe the basic ideas that we 
will use in our analysis. Although the same concepts are of use in three dimensional 
flows, we describe here just the two dimensional situation, since as discussed below the 
structures we identify can be considered two dimensional to a good approximation dur- 
ing the time scales relevant here. 

Stagnation points are well known features of steady flows that generally play an im- 
portant role in "organizing" the qualitatively distinct streamlines in the flow. For exam- 
ple, saddle-type stagnation points can occur on boundaries where streams of flow tan- 
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gential to the boundary coming from opposite directions meet and then separate from 
that boundary. Saddle-type stagnation points can occur in the interior of a flow at a point 
where fluid seems to both converge to the point along two opposite directions and di- 
verge from it along two different directions. In steady flows the stagnation point is a 
trivial example of a fluid particle trajectory (i.e., it is a "solution," in this case a fixed point 
solution, of the equations for fluid particle motions defined by the velocity field) and the 
saddle point nature is manifested by the fact that there are directions for which nearby 
trajectories approach the stagnation point at an exponential rate and move away from 
the stagnation point at an exponential rate. These directions are sometimes referred to as 
"stagnation streamlines". They define material curves that "cross" at the stagnation point 
and typically form the boundaries between qualitatively distinct regions of flow. 

A related "time-dependent" picture as that described above exists for unsteady flows, 
with both similar, as well as much more complex, implications for transport. In unsteady 
flows a stagnation point at a given time -or rather, an Instantaneous Stagnation Point 
(ISP)- is a location at which velocity vanishes at that time. The sequence of locations is 
generally not a fluid particle trajectory \lde et al.ll2002ft . The true analog of the sad dle-type 



stagn ation point of steady flows is a Distinguished Hyperbolic Trajectory (DHT) (|lde et al. 



2002). "Hyperbolic" is the dynamical systems terminology for "saddle-type". These are 



fluid particle trajectories that have (time-dependent) directions for which nearby trajec- 
tories approach and move away from the D HT a t expo nential rates. "Distinguished" 



is a notion that is discussed in detail in 



Ide et al 



(2002), but the idea is that these are 



the key, isolated hyperbolic trajectories that serve to organize the transport behavior in 
a flow because they remain substantially more localize d (in a well de f ined sense, se e 



Ide et al.l Il2002h') than n eighboring hyperbolic trajectories. Ilde et al . (2002). ITu et al . (2003), 



and 



Mancho et al. 



2004) develop the algorithms that allow us to compute DHTs in a given 
flow. They are iterative methods that start with a first guess for the DHT positions in an 
interval of time (i.e. an initial curve in space and time) and then refine the space-time 
curve by imposing the criteria of hyperbolicity and localization. In the flows considered 
here, a good first guess is the location of ISPs, since it turns out that a DHT is often found 
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in the neighborhood of an Eulerian ISP. 

Just as in the steady case, there are analogs to the stagnation streamlines: in the dy- 
namical systems terminology these are referred to as the stable and unstable manifolds of the 
DHT, and they are time-dependent material curves. In dynamical systems terminology 
the fact that they are material curves means that they are invariant, i.e. a fluid particle 
trajectory starting on one of these curves must remain on that material curve during the 
course of its time evolution. "Stable manifold" means that trajectories starting on this 
material curve approach the DHT at an exponential rate as time goes to infinity, and "un- 
stable manifold" means that trajectories starting on th is material curve approach the DHT 



Mancho et al. 



develop 



at an exponential rate as time goes to minus infinity, 
the algorithms that enable us to compute the stable and unstable manifolds of hyperbolic 
trajectories. 

In unsteady flows stable and unstable manifolds of DHTs can intersect in isolated 
points different from the DHTs. This is a fundamental difference with respect to steady 
flows, and give rise to moving regions of fluid bounded by pieces of stable and unstable 
manifolds, the so called "lobes". Since the manifolds are material lines, fluid can not cross 
them by purely advective processes and thus they are perfect Lagrangian barriers (diffu- 
sion, or motion along the third dimension can however induce cross-manifold transport). 
Motion of the "lobes" is thus the mechanism responsible for mediating Lagrangian trans- 
port between different regions . References describing "lobe dynamics" in general are 



Rom- 



1998); 



Kedar and Wiggins ( 199(1 feeigie etal 1994 ): lwigginslil992|) : Malhotra and Wiggins 



Samelson a nd Wiggins (2006). Exam 
Ngan and Shepherd (199 



graphic flows are 



Miller et al 



(2002); 



es of applications of l o be dynamics to oceano 



Rogerson et al.l(|199' 



1 



Yuan et al 



(200llEo04); 



Deese et alJ (|2002|) . We will describe these ideas more fully in the con- 



text of transport associated with the Balearic front in the Mediterranean. 
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3 The Ocean Circulation Model 



In this work we analyze velocity fields which are obt ained from an ocean model, DieC AST 



1997|) , adapted to the Mediterranean Sea pietrich et al.l |2004; 



Fernandez et al. 



Dietrich 

2005^ . The 3D primitive equations are discretized with a fourth-order collocated control 



volume method. In zones adjacent to boundaries a conventional second order method is 
used. A fundamental feature of control volume based models is that the predicted quan- 
tities are control volume averages, while face-averaged quanti t ies ar e used to evaluate 
fluxes across control volume faces (panderson and Brassingtonl |1998). These quantities 
are computed using fourth-order approximations and numer ical dispersion errors are 
further reduced in the modified incompressibility algorithm by Di etrichl (|1997|) . 

Horizontal resolution is the same in both the longitudinal (4>) and latitudinal (A) direc- 
tions, with A$)= (1/8) of degree and AA = cos AA(J) thus making square horizontal control 
volume boundaries. Vertical resolution is variable, with 30 control volume layers. The 
thickness of control volumes in the top layer is 10.3 m and they are smoothly increased 
up to the deepest bottom control volume face at 2750 m. Thus ETOP05 bathymetry is 
truncated at 2750 m depth and it is not filtered or smoothed. 

Horizontal viscosity and diffusivity values are constant and equal to 10 m 2 s _1 . For 
the vertical viscosity and diffusivity, a formu lation based on the Richardson number de- 
veloped in Pacanowski and Philander jl981 ) is used, with background values set at near- 
molecular values (1 0~ 6 and 2 x 1 0~ 7 m 2 s _1 respectively). We use monthly mean wind 
stress reanalyzed from 10 m wind ou tput from ECMWF, a s chosen for the Mediterranean 



Sea Models Evaluation Experiment (peckers et al. 



2002). The heat and the freshwater 



fluxes used to force the model are model-de termined from monthly climatological SST 
and SSS as described in 



Dietrich et al 



(2004). The only open boundary is the Strait of 
Gibraltar, where inflow conditions are set similar to observations and outflow is model- 
determined by upwind. Everywhere else, free-slip lateral boundary conditions are used. 
All bottom dissipation is represented by a conventional nonlinear bottom drag with a co- 
efficient of 0.002. Lateral and bottom boundaries are thermally insulating. The model is 
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initialized at a state of rest with the annual mean temperature and salinity fields taken 
from the climatological data. The spin-up phase of integration is carried out for 16 years. 
Each year is considered to have 12 months 30 days length each (i.e. 360 days). The cli- 
matological forcings we use are adequate to identify the mechanisms and processes oc- 
curring under typical or average circumstances. Under this approach, high frequency 
motions are weak in our model and a daily sampling is adequate. The impact on trans- 
port of disturbances containing high frequencies, such as storms or wind bursts, is not 
the focus of the present Paper and would need specific modelling beyond climatological 
forcing. 

We focus on velocity fields obtained at the second layer which has its center at a depth 
of 15.93 m. This is representative of the surface circulation and is not as directly driven 
by wind as the top layer. We have recorded velocities, temperatures and salinities in 
this model layer for five years. Dynamical systems approaches have already been ap- 
plied to this data set , in particular Lyapunov techn iques to quantify mix ing strength 
( d'Ovidio et al. 12004 ), and the "leaking" approach ( Schneider et al. 2005 ) to quantify 
escape and residence times in several areas. Here we concentrate in the Northwestern 
region and apply the methods of lobe dynamics to characterize transport processes in the 
North Balearic Front area. 

Figure [l] shows an example of the output of the model for the velocity field in the 
selected layer of the Western Mediterranean Sea at day 649 (the 19th day of the tenth 
month -October- of the second year). Two well known currents, the Northern Current 
flowing southwards close to the Spanish coast and the Balearic Current associated with 
the North Balearic Front and flowing northeastwards North of the Balearics, are observed 
although significantly deformed by the presence of eddies. The Figure shows also the ISPs 
of the velocity field: circles for the elliptic and crosses for the saddle-type ones. 

The velocity field has small vertical components, so that this is not strictly a two di- 
mensional flow. With vertical velocities of the order of 1 5 m/ s, particles in the sec- 
ond model layer require about 13 days to traverse the layer. But during that time this 
vertical velocity is not constant. Averaging over the relevant time scales we find effec- 
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Figure 1: Velocity field at 15.93 m depth in the Northwestern Mediterranean at simulation 
day 649 (October). The Westernmost coast is the Spanish one, the islands are the Balearics, 
and the Easternmost coasts are portions of Corsica and Sardinia. Circles indicate elliptic 
ISPs and crosses indicate saddle-type ISPs. 
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tive velocities o f 0.1-0 .7 m/ day, depending on location and season (jd'Ovidio et al.l 12004 : 



Schneider et al. 



200*5!) . an d thus residence times in the second layer are between two 



weeks and several months. As a rule of thumb we can consider that trajectories preserve 
two dimensionality during time intervals of about 20 days. Since most of our trajectory 
integrations will be restricted to time intervals below that duration, they can be consid- 
ered two dimensional to a good approximation. 

By inspection of the temperature and salinity model outputs we identify a four month 
interval starting on October of the second simulation year (i.e. Autumn and early Winter) 
as an example in which gradients are strong most of the time and thus the North Balearic 
Front well defined (see Fig. EJ. We select this interval of time as our study case for which 
dynamic structures will be calculated. 



4 Computation of Trajectories and Manifolds in Mediter- 
ranean Data Sets 

The equations of motion that describe the horizontal evolution of particle trajectories in 
our velocity field are 

u(4»,A,t) 



d4> 
dt 
dA 
dt 



Rcos(A) 

yjjvVtj 

R 



(i) 

(2) 



where u and v represent the eastwards and northwards components of the surface veloc- 
ity field coming from the simulations described in the previous section, R is the radius of 
the Earth (6400 km in our computations), <J) is longitude and A latitude. Particle trajec- 
tories must be integrated in equations (l)-(2) and since information is provided just in a 
discrete space-time gri d, a first issue to dea l with is that of interpolation of discrete data 
sets. A recent paper by iMancho et alJ (|2006|) compares different interp olation techniques 
in tracking particle trajectories. Bicubic spatial interpolation in space (|Press et al.l Il992|) 
and third order Lagrange polynomials in time are shown to provide a computationally 
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efficient and accurate method. We use this technique in our c omputations. How ever we 



Press et al. 



1992) requires 



notice that bicubic spatial interpolation in space as discussed in ] 
an equally spaced grid. Our data input is expressed in spherical coordinates, and the grid 
is not uniformly spaced in the latitude coordinate. In order to interpolate in an uniformly 
spaced grid, we transform our coordinate system (A, <£) to a new coordinate system with 
coordinates 40- The latitude A is related to the new coordinate \x 



\x = ln|secA + tanA| 



(3) 



Our velocity field is now on a uniform grid in the 4>) coordinates. The equations of 
motion in the old variables are transformed to a new expression in the new variables, 

u(4>, n,t) 



deb 
dt 
chx 
dt 



Rcos(A(n)) 



Rcos(A(n)) 

where A(|x) is obtained by inverting Eq. ©, i.e. 

A 



7T 



- -2atan(e~^) 



(4) 
(5) 

(6) 



Once trajectories are integrated from these equations, for presentation purposes one 
can convert \x values back to latitudes A just by using ©. 

Distinguished hyperbolic trajectories and their unstable and stable manifolds are the 
main dynamical systems objects that we use to describe and quantify transport. Our ve- 
locity data set lasts only for a finite time, it is highly aperiodic in time, and turbulent - 
Comput ation of hyper bolic trajectories for this kind of flow is discussed in lTu et al. 



and in 



Ide et al 



2002), where algorithms for their computation are developed. A novel 



technique to compute sta ble and unstab l e man ifol ds of hyperbo l ic tr ajectories in aperi 
odic flows is developed in 



Mancho et al. 



fl2003D . In lMancho et all (12004)) these hyperbolic 
trajectory and manifold computation algorithms are combined into a unified algorithm 
and successfully applied to a turbulent wind driven, quasigeostrophic double-gyre. We 
will now apply these same algorithms to compute hyperbolic trajectories and their stable 
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and unstable manifolds that we will use to describe and quantify transport associated 
with the North Balearic front. 



5 Lagrangian structures and transport in the Balearic Sea 



We focus on the region North of the Balearic Islands, the Balearic Sea. The main oceanic 
structures kno wn to be present there are the Balearic current and the associated North 



Balearic Front (Lo pez-Garcia et al.l [1994; 



Millot 



1999). This last feature is known to rep- 



resent a transition zone between saltier and fresher waters in the Western Mediterranean. 
The salinity fields obtained from our computer simulation display significant salinity gra- 
dients (and also temperature gradients in Winter-in Summer the surface layer is heated 
in a rather homogeneous way) in the area (see Fig. |2j|. Our aim here is to interpret the 
presence of the gradients and the front in terms of a semipermanent "Lagrangian barrier" 
across which little transport occurs. This construction would also reveal the routes along 
which this transport happens. Topological changes in that picture, associated with the 
crossing by eddies and that may be interpreted as the breakdown of the front, are also 
observed during the simulation. 

Gradients are rather well defined during Autumn of the second simulation year and 
early Winter of the third one. During this period we find a long interval (from day 
649 to day 731) in which a Lagrangian structure constructed using stable and unstable 
manifolds of DHTs remains persistent and acts as a partial barrier to transport. Its lo- 
cation is well correlated with the salinity gradients, so that it can be interpreted as a 
Lagrangian identification of the North Balearic Front. The weak transport across the 
st ructure can be described in te rms of lobe dynamics. The situation resembles the one 



m 



Coulliette and Wiggins (2001) in which quasigeostrophic dynamics has been used to 



model a double gyre situation and the central jet between the two gyres in that problem 
plays a similar role as that of the Balearic current in our problem. Lobe dynamics was 
there successfully applied to quantify transport across the jet, occurring by the so called 
turnstile mechanism. However in the more realistic data set analyzed here we need to 
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generalize some ideas used there. For example, in ICoulliette and Wiggins! ([2001) DHTs 
stay on one dimensional boundaries for all time (as boundaries are invariant), however in 
our scenario the relevant current does not start nor end on clear one dimensional bound- 
aries. This introduces some ambiguity in the identification of the relevant DHTs (and of 



the s addle-type ISPs used as starting positions in the iterative algorithm of l|Mancho et al 



2004) that we use to determine the DHTs) from which to compute the manifolds that will 
define the Lagrangian barrier. Most of the time, however, pieces of manifolds computed 
from different DHTs in the same area rapidly converge towards each other, thus indicat- 
ing that the location of the dominant hyperbolic curve is not a property of the particular 
choice of DHTs, but a property of the flow. Changes in the topology of the flow cause this 
convergence property to be lost. This happens at the end of the interval of time chosen in 
our study case and will be commented on below. 



5.1 Cross frontal transport: the turnstile mechanism 

The geometric objects used to chara cterize transport by lobe dynamics methods are con 
structed by the rules discussed in Malhotra and Wiggins! 1998 ; Coulliette and Wiggins 



2001 



Samelson and W iggins 2006). Here we describe in some detail these procedures in 



a way that is particular to our flow situation. Figs. |2]-|4| illustrate the construction at days 
649 and 657. Fig. |3 contains all the dynamic structures superimposed on a salinity field, 
whereas for clarity only the 'boundary' and the 'lobes' are displayed in Figs. |3j and HJ 
respectively. 

1. Two DHTs should be identified, one in the western part of the front to be character- 
ized and another in the eastern, that persists close to their initial positions during 
the whole time interval of interest (denoted by [to, tisi]). Since our algorithm to locate 
DHTs uses saddle ISP positions as first guesses, figures displaying ISPs such as Fig. 
|T]are used to estimate these positions and the temporal persistence of the ISPs. The 
positions calculated for the selected DHTs are plotted in Figs. |2]-|6]as black dots and 
labelled as H w (the western one) and H E (the eastern one). 
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2. As the mean current flows eastwards we proceed as in Coulli ette and Wi ggins ( 2001) 
and compute the unstable manifold of the western DHT and the stable manifold of 
the one in the East. They are the red and blue lines in Fig. |2j respectively. For clarity 
in the presentation, of the two branches of each manifold (one at each side of the 
DHT from which they emanate) we display in Fig. |2] only the one pointing in the 
direction of the other DHT. Along both pieces of manifolds, in the region between 
the two DHTs, the dominant direction of motion is from west to east. As character- 
istic in unsteady flows, both manifolds intersect repeatedly. Some of the intersection 
points are marked with cyan dots. It is also typical that the unstable manifold (red) 
of a DHT (Hw in this case) emerges from it relatively straight, whereas it fluctuates 
widely when approaching the vicinity of the opposite DHT. In the same way, the 
stable manifold (blue) of H E joins it smoothly, but displays characteristic oscillations 
when close to the western DHT. The stable and unstable manifolds displayed in Fig. 
12 a) have been computed after backwards and forwards time integration of a small 
segment in the direction of the stable and unstable linear subspaces in the neigh- 
borhood of the DHT for time periods of 14 and 19 days, respectively. Integrations 
for longer time periods provide longer manifolds. However due to the restrictions 
of the two dimensional approximation, a time period beyond 20 days is not com- 
pletely trustworthy. In practice this means that the pieces of the manifold that are 
far from the DHT may deviate from the true manifold as those pieces are the ones 
obtained with longer time integrations. For instance, the unstable manifold in Fig. 
|2]b) has been computed for a time integration of 27 days. This means that unstable 
structure in the far East is not completely reliable. However the piece of manifold 
governing the turnstile mechanism, which is closer to the DHT, is obtained with 
numerical integrations below the validity limit of 20 days and predictions obtained 
from them are correctly approached in the two dimensional approximation. 

3. Between the beginning and the end of the chosen time interval, we choose a se- 
quence of times at which to analyze the manifold positions and compute objects 
relevant for Lagrangian transport. The sequence of "observation times" is denoted 
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by t < t] < t 2 < . . . < t N _T < t N . We note that the ti do not need to be equally 
spaced. In order to illustrate the construction of boundaries and the turnstile mech- 
anism for crossing the boundaries we need only two times, and for this purpose we 
will choose to show days 649 and 657. 

4. At each of the selected times t t , a "boundary" is constructed by choosing a finite 
piece of the unstable manifold of the western DHT, [H w , b t . ], and a finite piece of the 
stable manifold of the eastern DHT, [b ti , H E ], so that they intersect in precisely one 
point, b ti , which is called a boundary intersection point. The points {b t J, in addition to 
satisfying an ordering constraint specified below, should be selected in such a way 
that a boundary relatively straight, i.e. free of the violent oscillations displayed by 
each of the manifolds when approaching the opposite DHT, is obtained. Since the 
boundary is pinned at the points H w and H E we obtain a sequence of boundaries 
that fluctuate in time but remain approximately in the same place. Figure |2] shows 
the selection of the boundary intersection points b 649 and b 657 at times 649 and 657, 
respectively. For clarity, Fig. |3]displays just the resulting boundaries. 

Since the boundary is made of material lines, no fluid can cross it by horizontal ad- 
vection processes, except at the observation times ti , t 2 , . . . , t N at which the boundary 
is redefined. At these times, the only way fluid at one side of the boundary can be 
transferred to the other side is by the turnstile mechanism described and quanti- 
fied below. If this transport amount is small (as it will be shown to be the case), 
the boundary can be characterized as a "barrier", and gradients will be maintained 
across it. As seen in Fig. |2J the position of the boundary is well correlated with the 
position of the salinity front thus confirming that the dynamical systems techniques 
developed here are useful to identify the North Balearic front in terms of Lagrangian 
objects. 

5. Construct turnstiles at ti. At time tt we consider the point, denoted by b^7 +1 , that 
will evolve into the boundary intersection point b ti+1 (clearly, we cannot do this at 
t N ). Since the stable and unstable manifolds are invariant, this point is also on both 
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the stable and unstable manifolds. In the same way, bi is the location at time tt of 
the boundary intersection point that was located at b tt _, at the previous time t t n . 
The additional constraint that needs to be imposed when choosing the sequence 
{b t J is that b^ 1 results "upstream" (i.e. closer to H w along its unstable manifold, 
or further from H E along its stable manifold) from b ti . This introduces a restriction 
on the choice of b tl+1 once b ti is chosen. Since ordering of points along manifolds 
is preserved by time evolution, it turns out that bi would be "downstream" from 
b t .. Figures El and H| display some of these intersection points, showing that our 
selection satisfies the ordering constraints. The segments of stable and unstable 
manifolds between b ti and b^ +1 (and between b ti and b£_.,) trap regions of fluid, 
and these regions of fluid defined in this way are referred to as lobes. Some of them 
can be seen in Fig. |2j and more clearly in Fig. |U 

6. Consider the evolution of the turnstile lobes from ti to t| + i . In this case b^ +1 evolves 
to the boundary intersection point, b tt+1 , and the boundary intersection point at tt 
evolves to a point b£, which is also on both the stable and unstable manifolds. Then 
the lobes between b ti+1 and b^ represent the time evolution of the turnstile lobes 
from ti to tj, + i . Note that, because of the way in which the boundary is redefined 
at each observation time, and because of the different shape of the manifolds when 
close or far from the DHT from which they emanate, each lobe is on opposite sides 
of the boundary at the two considered times. These turnstile lobes contain all the 
fluid that has crossed the boundary between ti and t i+1 . This transport process is 
illustrated more clearly in Fig. |4] where the lobes experiencing the turnstile mecha- 
nism are plotted before and after crossing the boundary. 

The geometric construction performed at every observation time tt as explained above 
allows us to calculate the amount of transport occurring across the boundary during each 
time interval. 



17 



In computing the area A(F) of a lobe V we use the formula 



A(D 



sinAdcj) 



(7) 



ar 



where the integration is around the closed curve which forms the boundary of the lobe. 
That this formula gives the area can b e seen by co nsidering the differential form tu = 
— R 2 siuAd(j), calculating its differential ( Spivakl 119650 : do; = R 2 cosAdcj)dA which is iden- 
tical to the area ele ment on the sphere in spherical coordinates, and recalling Stokes theo- 



rem 



jspivak 



1965): 



ar 



duo. 



(8) 



For day 649, as shown in Fig.HJ the two turnstile lobes have areas of 493.2km 2 (the lobe 
below the boundary, to the east) and 71 6.9km 2 (the lobe to the west above the boundary). 
At day 657, the eastern lobe is above the new boundary, and the western lobe is below it. 
Assuming that the divergence of the surface flow can be neglected so that the areas are 
unchanged - in numerical experiments we have never observed more than a 3% change 
- we can calculate the flux across the barrier to be (71 6.9 — 493.2)km 2 = 223.7km 2 , in the 
southerly direction. 

Modified Atlantic Waters occupy the surface layers of the area until an average depth 
of about 1 50m. Using as an approximation for the mean horizontal speed of this water 
mass the values at the second model layer considered here, multiplication of the 1 50m 
depth by the area of the lobes gives 33.56 x 1 9 m 3 in 8 days, or an average flux over this 
interval of 0,049 Sv (lSv=10 6 m 3 /s). The average flux obtained from area calculations of 
further turnstile lobes until the middle of November remains below that value (the aver- 
age is 0.025 Sv, always southwards). This should be compared to the 0.75-0.5 Sv which are 
transported by the Balearic current. We see that cross-boundary transport is small dur- 
ing this time interval and thus the Lagrangian boundary acts as a "barrier to transport" 
permitting only small amounts of mixing between Northern and Southern waters. It will 
maintain a salinity (and thus density) front that we identify with the observed North 



18 



Balearic front. There is some indeterminacy in the definition of the boundary that we 
identify as the front, arising from some freedom in the selection of the the intersections 
{b t J. But other choices can only displace the boundary by a distance of the order of the 
size of the lobes, which we see is small when not too close to the DHTs, and in fact also of 
the order of the width of the transition region in salinity distributions such as the one in 
Fig. |2j i.e. the width of the front. 

5.2 Spatio-temporal structure of cross frontal transport 

Since the only way in which our constructed Lagrangian boundaries can be crossed is via 
the turnstile mechanism, the earlier and later location of the turnstile lobes reveals the 
dominant routes along which the weak cross-front transport occurs. Figure |5] shows the 
time evolution of a pair of turnstile lobes, one initially above and the other initially below 
the boundary, as they evolve in time. The crossing of the boundary by the turnstile mech- 
anism occurs between days 676 and 681, and the remaining panels in the figure show the 
position of these lobes at some earlier and later times. The sequence illustrate that lobes 
move essentially along the boundary except when close to H E , where they are ejected as 
filaments transverse to the front (the one initially in the south ejected towards the north 
and reciprocally for the one started in the north) and when close to H w , where they have 
also the shape of transverse filaments and they become entrained into the boundary re- 
gion. Fig. |5]also illustrates how lobes transport water of different salinity (coded in colors) 
and how the above routes for lobe motion and shape correlate with the salinity distribu- 
tion in the area. Note that the length of the manifolds and the whole process depicted in 
Fig. |5] lasts only 21 days, so that the plotted manifolds remain approximately horizontal, 
with only small corrections from the vertical flow. 

5.3 An Eddy-front interaction: disruption of the Lagrangian boundary 

Not all flow configurations allow the geometric construction identifying the Lagrangian 
boundary and the associated turnstile lobes to be performed. It may happen that no pair 
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of DHTs persist in a given area long enough to support the mechanism, or their manifolds 
can fail to intersect. It may also happen that manifolds started at relatively close DHTs do 
not converge to the same curve but remain significatively distinct, so that a unique well 
defined boundary can not be properly identified. In such cases, the turnstile transport 
mechanism is not the most relevant one. At the end of the simulation interval analyzed 
here we observe a change in the topology that signals the end of the predominance of the 
turnstile mechanism in a process that can be interpreted as the breakdown of the front by 
the interaction with an eddy. As a first symptom, calculations of turnstile lobe areas reveal 
an increase in cross-front transport starting at day 674 (middle of November). The average 
transport between days 674 and 700 is of 0.303 Sv (southwards), still smaller than the 
Balearic current transport but significantly larger than the average cross-frontal transport 
during the previous month (0.025 Sv, also southwards). At day 711, stable manifolds 
emerging from He and from another rather close DHT cease to converge into each other, 
signaling the end of a situation with an essentially unique well defined boundary. Later, 
at day 731 our algorithm is unable to find the location of H E starting from ISPs in the area. 
This probably means that H E has moved away from the area under study. The black dot 
in Fig. [6]is another DHT found in the region. But its stable manifold (i.e., the finite length 
of manifold that we are able to compute) does not intersect the unstable manifold from 
H w , thus revealing that is in fact a DHT different from H E , and that it can not support the 
turnstile mechanism. The time evolution of the unstable manifold from H w suggests that 
the reason from the change in behavior is the breakdown of the Lagrangian barrier by 
the crossing of an eddy, identified by the rolling up of the manifold around an elliptic ISP 
(Fig. |5J. Note that even in this situation the manifold position is well correlated with the 
salinity distribution, thus indicating that still Lagrangian structures are relevant. But the 
transport mechanism is clearly different from the turnstile described above, being more 
appropriately described as water transport inside an eddy. 
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6 Conclusions 



In this work we have applied in a systematic way some tools developed in the context 
of dynamical systems theory and known generically as 'lobe dynamics'. The computer 
generated surface velocity field studied here is more complex and less regular that other 
velocity fields previously considered in this context, but we have found that one of the 
main mechanisms of transport by lobe motion, the turnstile, is still at work. The method- 
ology includes the construction of a 'barrier' across which to compute transport, and in 
our application to the Northwestern Mediterranean dynamics it has been identified with 
one of the main oceanographic structures present here, the North Balearic Front. Trans- 
port across it proceeds in the form of filaments that are entrained into the front close to 
a DHT 'upstream', and released also in the form of filaments close to another DHT lo- 
cated 'downstream'. The ejection of these filaments at that location can ex plain recent 
observations of waters saltier than expected just east of the island of Menorca (Emelianov 



2006). The identification of the DHTs and the calculation of their locations is by itself an 



important subject, since they organize the flow in the area and, because of this and of the 
sensitivity of the trajectories in their neighborhood, they are candidates for launch loca 
tions in efficiently designed drifter release experiments (|Poje et al.l |2002; 



Molcard et al 



2006). 



Despite the success of the approach described here, much work remains to be done 
in order to develop dynamical systems techniques into a collection of systematic tools 
for analyzing general oceanographic data. A classification and understanding of the dif- 
ferent topological regimes leading to qualitatively different modes of transport and the 
transitions among them, similar to the existing ones for steady and time-periodic flows, 
would be desirable for the cases of turbulent aperiodic flows. A characteristic of the dy- 
namical systems approach is that it provides an unusually high detailed description of the 
spatio-temporal structure of Lagrangian transport. Therefore it may well turn out to be 
the optimal tool for analyzing data from ocean models with much higher spatio-temporal 
resolution (e.g., higher frequency atmospheric forcing and more resolved spatial scales) 
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which capture more physical processes. This would enable us to better define, for exam- 
ple, the process of the destruction of a barrier. In addition, consideration of the impact of 
non-Lagrangian processes such as diffusion, of vertical motions, and of strong localized 
perturbations beyond climatological forcing such as storms, would be needed to have a 
more complete vision of transport phenomena and mechanisms. 

Acknowledgements 

A.M.M. acknowledges the MCyT (Spanish Government) for a Ramon y Cajal Research 
Fellowship and financial support from MEC (Spanish Government) reference MTM2004- 
00797 and the Royal Society-CSIC cooperation agreement reference B2003GB03. E.H.- 
G. acknowledges financial support from MEC and FEDER through Project CONOCE2 
(FIS2004-00953). We also acknowledge M. Emelianov for communicating us results from 
a recent cruise before publication. D.S. and S. W. acknowledge financial support from 
ONR Grant No. N00014-01-1-0769. 

References 

Bartlett, J. H., 1982: Limits of stability for an area-perserving polynomial mapping. Cel. 
Mech., 28, 295-317. 

Beckers, J.-M., M. Rixen, P. Brasseur, J.-M. Brankart, A. Elmoussaoui, M. Crepon, Herbaut, 
F. Martel, F. V. den Berghe, L. Mortier, A. Lascaratos, P. Drakopoulos, G. Korres, K. Nit- 
tis, N. Pinardi, E. Masetti, S. Castellari, P. Carini, J. Tintore, A. Alvarez, S. Monserrat, 
D. Parrilla, R. Vautard, and S. Speich, 2002: Model intercomparison in the Mediter- 
ranean: MEDMEX simulations of the seasonal cycle. /. Mar. Sys., 33, 215-251. 

Beigie, D., A. Leonard, and S. Wiggins, 1994: Invariant manifold templates for chaotic 
advection. Chaos, Solitons, and Fractals, 4(6), 749-868. 



22 



Channon, S. R. and J. L. Lebowitz, 1980: Numerical experiments in stochasticity and ho- 
moclinic oscillations. Ann. New York Acad. Sci, 357, 108-118. 

Coulliette, C. and S. Wiggins, 2001: Intergyre transport in a wind-driven, quasi- 
geostrophic double gyre: An application of lobe dynamics. Nonlinear Processes in Geo- 
physics, 8, 69-94. 

Deese, H. E., L. J. Pratt, and K. R. Helfrich, 2002: A laboratory model of exchange and 
mixing between western boundary layers and subbasin recirculation gyres. /. Phys. 
Oceanogr., 32(6), 1870-1889. 

Dietrich, D., 1997: Application of a modified "a" grid ocean model having reduced nu- 
merical dispersion to the Gulf of Mexico circulation. Dyn. Atmos. Oceans, 27, 201-217. 

Dietrich, D., R. Haney, V. Fernandez, S. Josey, and J. Tintore, 2004: Air-sea fluxes based 
on observed annual cycle surface climatology and ocean model internal dynamics: a 
non-damping zero-phase-lag approach applied to the Mediterranean sea. /. Mar. Sys., 
52, 145-165. 

d'Ovidio, R, V. Fernandez, E. Hernandez-Garria, and C. Lopez, 2004: Mixing structures 
in the Mediterranean sea from finite-size Lyapunov exponents. Geophys. Res. Lett., 31, 
L12203 (1-4), doi:10.1029/2004GL020328. 

Duan, J. Q. and S. Wiggins, 1996: Fluid exchange across a meandering jet with quasi- 
periodic time variability. /. Phys. Oceanogr., 26, 1176-1188. 

Emelianov, M., 2006: Personal communication. 

Fernandez, V., D. Dietrich, R. Haney, and J. Tintore, 2005: Mesoscale, seasonal and inter- 
annual variability in the Mediterranean sea using a numerical ocean model. Progress in 
Oceanography, 66, 321-340. 



23 



Ide, K., D. Small, and S. Wiggins, 2002: Distinguished hyperbolic trajectories in time de- 
pendent fluid flows: analytical and computational approach for velocity fields defined 
as data sets. Nonlinear Processes in Geophysics, 9, 237-263. 

Isern-Fontanet, J., Garcia-Ladona, and J. Font, 2006: The vortices of the Mediterranean 
sea: an altimetric perspective. /. Phys. Ocean., 36, 87103. 

Jones, C. K. R. T. and S. Winkler, 2002: Invariant manifolds and Lagrangian dynamics in 
the ocean and atmosphere. Handbook of dynamical systems, North-Holland, Amsterdam, 
55-92. 

Ju, N., D. Small, and S. Wiggins, 2003: Existence and computation of hyperbolic trajecto- 
ries of aperiodically time-dependent vector fields and their approximations. Int. J. Bif. 
Chaos, 13, 1449-1457. 

Lopez-Garcia, M., C. Millot, J. Font, and E. Garcia-Ladona, 1994: Surface circulation vari- 
ability in the Balearic Basin. /. Geophys. Res., 99 (C2), 3285-3296. 

MacKay, R. S., J. D. Meiss, and I. C. Percival, 1984: Transport in Hamiltonian systems. 
Physica D, 13, 55-81. 

Malhotra, N. and S. Wiggins, 1998: Geometric structures, lobe dynamics, and Lagrangian 
transport in flows with aperiodic time-dependence, with applications to Rossby wave 
flow. /. Nonlinear Science, 8, 401-456. 

Mancho, A. M., D. Small, and S. Wiggins, 2004: Computation of hyperbolic and their sta- 
ble and unstable manifolds for oceanographic flows represented as data sets. Nonlinear 
Processes in Geophysics, 11, 17-33. 

— 2006: A comparison of methods for interpolating chaotic flows from discrete velocity 
data. Computers & Fluids, 35, 416-428. 



24 



Mancho, A. M., D. Small, S. Wiggins, and K. Ide, 2003: Computation of stable and 
unstable manifolds of hyperbolic trajectories in two-dimensional, aperiodically time- 
dependent vector fields. Physica D, 182, 188-222. 

Miller, P. D., L. J. Pratt, K. R. Helfrich, and C. K. R. T. Jones, 2002: Chaotic transport of 
mass and potential vorticity for an island recirculation. /. Phys. Oceanogr., 32(1), 80-102. 

Millot, C, 1999: Circulation in the western Mediterranean sea. /. Mar. Sys., 20, 423-442. 

Molcard, A., A. Poje, and T. Ozgokmen, 2006: Directed drifter launch strategies for La- 
grangian data assimilation using hyperbolic trajectories. Ocean Modell., 12, 268-289. 

Ngan, K. and T. G. Shepherd, 1997: Chaotic mixing and transport in Rossby wave critical 
layers. /. Fluid. Mech., 334, 315-351. 

Olson, D., 1991: Rings in the ocean. Annu. Rev. Earth Planet. Sci, 19, 283-311. 

Pacanowski, R. C. and S. G. H. Philander, 1981: Parametrization of vertical mixing in 
numerical models of tropical oceans. /. Phys. Oceanog., 11(11), 1443-1451. 

Pinot, J. M., J. Tintore, and D. Gomis, 1995: Multivariate analysis of the surface circulation 
in the Balearic sea. Prog. Oceanog., 36, 343-376. 

Poje, A. C, M. Toner, A. D. Kirwan, and C. K. R. T. Jones, 2002: Drifter launch strategies 
based on Lagrangian templates. /. Phys. Oceanogr., 32(6), 1855-1869. 

Press, W. H., S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, 1992: Numerical Recipes 
in C. Cambridge University Press. 

Puillat, I., I. Taupier-Letage, and C. Millot, 2002: Algerian eddies lifetime can be near 3 
years. /. Mar. Sys., 31, 245-259. 

Rogerson, A. M., P. D. Miller, L. J. Pratt, and C. K. R. T. Jones, 1999: Lagrangian motion 
and fluid exchange in a barotropic meandering jet. /. Phys. Oceanogr., 29, 2635-2655. 



25 



Rom-Kedar, V., A. Leonard, and S. Wiggins, 1990: An analytical study of transport, mix- 
ing, and chaos in an unsteady vortical flow. /. Fluid Mech., 214, 347-394. 

Rom-Kedar, V. and S. Wiggins, 1990: Transport in two-dimensional maps. Arch. Rat. Mech. 
Anal, 109, 239-298. 

Ruiz, S., J. Font, M. Emelianov, J. Isern-Fontanet, C. Millot, J. Salas, and I. Taupier-Letage, 
2002: Deep structure of an open sea eddy in the Algerian Basin. /. Mar. Sys., 33-34, 
179-195. 

Samelson, R. and S. Wiggins, 2006: Lagrangian Transport in Geophysical Jets and Waves. 
Springer- Verlag, New York. 

Samelson, R. M., 1992: Fluid exchange across a meandering jet. /. Phys. Oceanogr., 22, 
431-440. 

Sanderson, B. and G. Brassington, 1998: Accuracy in the context of a control-volume 
model. Atmosphere-Ocean, 36, 355-384. 

Schneider, J., V. Fernandez, and Hernandez-Garcia, 2005: Leaking method approach to 
surface transport in the Mediterranean sea from a numerical ocean model. /. Mar. Sys., 
57, 111-126. 

Spivak, M., 1965: Calculus on Manifolds. Perseus Books, New York. 

Wiggins, S., 1992: Chaotic Transport in Dynamical Systems. Springer- Verlag, New York. 

— 2005: The dynamical systems approach to Lagrangian transport in oceanic flows. Annu. 
Rev. Fluid Mech., 37, 295-328. 

Yuan, G.-C, L. J. Pratt, and C. K. R. T. Jones, 2001: Barrier destruction and Lagrangian 
predictabilty at depth in a meandering jet. Dyn. Atmos. Oceans, 35, 41-61. 

— 2004: Cross-jet Lagrangian transport and mixing in a 2\ layer model. /. Phys. Oceanogr., 
34, 1991-2005. 



26 




IE 2E 3E 4E 5E 6E 7E 




IE 2E 3E 4E 5E 6E 7E 

Ion 



Figure 2: Salinity front and manifolds at days 649 and 657. Salinity (in psu) is color coded 
as indicated by the color bar. The black circles denote the DHT in the West (H w ) and 
the DHT in the East (H E ). These are the DHTs whose unstable (red) and stable (blue) 
manifolds, respectively, are used to construct the Lagrangian boundary of the front. The 
blue dots are boundary intersection points, as described in the text. 
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Figure 3: Boundaries at days 649 and 657 constructed from a (finite length) segment of 
the unstable manifold of H w and a (finite length) segment of the stable manifold of H E . 
The boundary intersection points are denoted by and b 657 , respectively 
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Figure 4: Turnstile lobes from day 649. There is precisely one intersection point between 
b649 and b^™ which implies that there are two lobes in the turnstile. The bottom figure 
shows the evolved position of these lobes at day 657. Comparing the two figures, one can 
see that the turnstile lobe to the North of the boundary at day 649 has moved to the South 
of the boundary at day 657, and the turnstile lobe to the South of the boundary at day 649 
has moved to the North of the boundary at day 657 (the location of the boundary itself 
can be seen in Fig. |3j. 29 



(a) Evolution of the turnstile lobes at day 674 
backwards in time to day 667, or, in other 
words, the location on day 667 of the lobes that 
are turnstile lobes at day 674. The notation 
bgg, denotes the evolution of the boundary in- 
tersection point at day 681 backwards two ob- 
servation times (to day 674 and to the current 
time, day 667). 
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(c) The turnstile lobes evolved to day 681. 
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(d) The turnstile lobes evolved to day 688. 
The notation bg^ denotes the evolution of the 
boundary intersection point at day 67 '4 forwards 
two observation times (to day 681 and to the 
current time, day 688). 



(b) Turnstile lobes at day 674. 



Figure 5: The spatio-temporal structure of cross frontal transport described by lobe dy- 
namics. Salinity (in psu) is color coded as indicated by the color bar. 
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(a) Day 733. The unstable manifold of Hw does 
not appear to intersect the stable manifold of 
H^. The unstable manifold is strongly influ- 
enced by the circulation associated with the el- 
liptic ISP immediately to the lower right of H^. 
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(c) Day 742. The earlier Lagrangian barrier 
now appears completely destroyed as a result 
of the influence of the eddy. 



(b) Day 736. As the elliptic ISP moves further 
away from its influence carries the unstable 
manifold with it. 

Figure 6: The mechanism of the disruption of the Lagrangian barrier by an eddy. Salinity 
(in psu) is coded in colors as indicated by the color bar. Red dots are ISPs of elliptic 
character. Green dots are saddle ISPs. Also shown as black dots are the DHTs Hw and 
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